wradlib logo png

wradlib data quality


Overview

Within this notebook, we will cover:

  1. Reading radar volume data into xarray based RadarVolume

  2. Wrapping numpy-based functions to work with Xarray

  3. Clutter detection

  4. Beam Blockage calculation

Prerequisites

Concepts

Importance

Notes

Xarray Basics

Helpful

Basic Dataset/DataArray

Matplotlib Basics

Helpful

Basic Plotting

Intro to Cartopy

Helpful

Projections

  • Time to learn: 10 minutes


Imports

import glob
import os

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import wradlib as wrl

Import Swiss Radar Data from CfRadial1 Volumes

We use some of the pyrad example data here. Sweeps are provided as single files, so we open each file separately and create the RadarVolume from the open Datasets.

fglob = "./data/example_pyrad/22179/MLL22179/MLL2217907250U*.nc"
flist = glob.glob(fglob)
flist.sort()
print("Files available: {}".format(len(flist)))
Files available: 20
ds = [xr.open_dataset(f, group="sweep_0", engine="cfradial1", chunks={}) for f in flist]
vol = wrl.io.RadarVolume(engine="cfradial1")
vol.extend(ds)
vol.sort(key=lambda x: x.time.min())
---------------------------------------------------------------------------
UFuncTypeError                            Traceback (most recent call last)
Cell In[5], line 1
----> 1 vol.sort(key=lambda x: x.time.min())

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/io/xarray.py:593, in XRadBase.sort(self, **kwargs)
    592 def sort(self, **kwargs):
--> 593     self._seq.sort(**kwargs)

Cell In[5], line 1, in <lambda>(x)
----> 1 vol.sort(key=lambda x: x.time.min())

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/_aggregations.py:1580, in DataArrayAggregations.min(self, dim, skipna, keep_attrs, **kwargs)
   1509 def min(
   1510     self,
   1511     dim: Dims = None,
   (...)
   1515     **kwargs: Any,
   1516 ) -> DataArray:
   1517     """
   1518     Reduce this DataArray's data by applying ``min`` along some dimension(s).
   1519 
   (...)
   1578     array(nan)
   1579     """
-> 1580     return self.reduce(
   1581         duck_array_ops.min,
   1582         dim=dim,
   1583         skipna=skipna,
   1584         keep_attrs=keep_attrs,
   1585         **kwargs,
   1586     )

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/dataarray.py:3681, in DataArray.reduce(self, func, dim, axis, keep_attrs, keepdims, **kwargs)
   3637 def reduce(
   3638     self: T_DataArray,
   3639     func: Callable[..., Any],
   (...)
   3645     **kwargs: Any,
   3646 ) -> T_DataArray:
   3647     """Reduce this array by applying `func` along some dimension(s).
   3648 
   3649     Parameters
   (...)
   3678         summarized data and the indicated dimension(s) removed.
   3679     """
-> 3681     var = self.variable.reduce(func, dim, axis, keep_attrs, keepdims, **kwargs)
   3682     return self._replace_maybe_drop_dims(var)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/variable.py:2014, in Variable.reduce(self, func, dim, axis, keep_attrs, keepdims, **kwargs)
   2012         data = func(self.data, axis=axis, **kwargs)
   2013     else:
-> 2014         data = func(self.data, **kwargs)
   2016 if getattr(data, "shape", ()) == self.shape:
   2017     dims = self.dims

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/duck_array_ops.py:382, in _create_nan_agg_method.<locals>.f(values, axis, skipna, **kwargs)
    380     with warnings.catch_warnings():
    381         warnings.filterwarnings("ignore", "All-NaN slice encountered")
--> 382         return func(values, axis=axis, **kwargs)
    383 except AttributeError:
    384     if not is_duck_dask_array(values):

File <__array_function__ internals>:200, in amin(*args, **kwargs)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/core.py:1748, in Array.__array_function__(self, func, types, args, kwargs)
   1745         return handle_nonmatching_names(func, args, kwargs)
   1747 if not hasattr(module, func.__name__):
-> 1748     return handle_nonmatching_names(func, args, kwargs)
   1750 da_func = getattr(module, func.__name__)
   1751 if da_func is func:

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/core.py:1727, in Array.__array_function__.<locals>.handle_nonmatching_names(func, args, kwargs)
   1724     args, kwargs = compute(args, kwargs)
   1725     return func(*args, **kwargs)
-> 1727 return _HANDLED_FUNCTIONS[func](*args, **kwargs)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/reductions.py:431, in min(a, axis, keepdims, split_every, out)
    428 @implements(np.min, np.amin)
    429 @derived_from(np)
    430 def min(a, axis=None, keepdims=False, split_every=None, out=None):
--> 431     return reduction(
    432         a,
    433         chunk_min,
    434         chunk.min,
    435         combine=chunk_min,
    436         axis=axis,
    437         keepdims=keepdims,
    438         dtype=a.dtype,
    439         split_every=split_every,
    440         out=out,
    441     )

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/reductions.py:235, in reduction(x, chunk, aggregate, axis, keepdims, dtype, split_every, combine, name, out, concatenate, output_size, meta, weights)
    232 else:
    233     reduced_meta = None
--> 235 result = _tree_reduce(
    236     tmp,
    237     aggregate,
    238     axis,
    239     keepdims,
    240     dtype,
    241     split_every,
    242     combine,
    243     name=name,
    244     concatenate=concatenate,
    245     reduced_meta=reduced_meta,
    246 )
    247 if keepdims and output_size != 1:
    248     result._chunks = tuple(
    249         (output_size,) if i in axis else c for i, c in enumerate(tmp.chunks)
    250     )

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/reductions.py:303, in _tree_reduce(x, aggregate, axis, keepdims, dtype, split_every, combine, name, concatenate, reduced_meta)
    301 if concatenate:
    302     func = compose(func, partial(_concatenate2, axes=sorted(axis)))
--> 303 return partial_reduce(
    304     func,
    305     x,
    306     split_every,
    307     keepdims=keepdims,
    308     dtype=dtype,
    309     name=(name or funcname(aggregate)) + "-aggregate",
    310     reduced_meta=reduced_meta,
    311 )

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/reductions.py:381, in partial_reduce(func, x, split_every, keepdims, dtype, name, reduced_meta)
    379 if is_arraylike(meta) and meta.ndim != len(out_chunks):
    380     if len(out_chunks) == 0:
--> 381         meta = meta.sum()
    382     else:
    383         meta = meta.reshape((0,) * len(out_chunks))

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/core/_methods.py:49, in _sum(a, axis, dtype, out, keepdims, initial, where)
     47 def _sum(a, axis=None, dtype=None, out=None, keepdims=False,
     48          initial=_NoValue, where=True):
---> 49     return umr_sum(a, axis, dtype, out, keepdims, initial, where)

UFuncTypeError: ufunc 'add' cannot use operands with types dtype('<M8[ns]') and dtype('<M8[ns]')
display(vol)
<wradlib.RadarVolume>
Dimension(s): (sweep: 20)
Elevation(s): (-0.2, 0.4, 1.0, 1.6, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 11.0, 13.0, 16.0, 20.0, 25.0, 30.0, 35.0, 40.0)
vol.root()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[7], line 1
----> 1 vol.root()

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/io/xarray.py:627, in RadarVolume.root(self)
    625 """Return root object."""
    626 if self._root is None:
--> 627     self.assign_root()
    628 return self._root

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/io/xarray.py:656, in RadarVolume.assign_root(self)
    652     sweep_fixed_angles = [ts.attrs["fixed_angle"].min() for ts in self]
    654 # extract time coverage
    655 times = np.array(
--> 656     [[ts.rtime.values.min(), ts.rtime.values.max()] for ts in self]
    657 ).flatten()
    658 time_coverage_start = min(times)
    659 time_coverage_end = max(times)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/io/xarray.py:656, in <listcomp>(.0)
    652     sweep_fixed_angles = [ts.attrs["fixed_angle"].min() for ts in self]
    654 # extract time coverage
    655 times = np.array(
--> 656     [[ts.rtime.values.min(), ts.rtime.values.max()] for ts in self]
    657 ).flatten()
    658 time_coverage_start = min(times)
    659 time_coverage_end = max(times)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/common.py:278, in AttrAccessMixin.__getattr__(self, name)
    276         with suppress(KeyError):
    277             return source[name]
--> 278 raise AttributeError(
    279     f"{type(self).__name__!r} object has no attribute {name!r}"
    280 )

AttributeError: 'Dataset' object has no attribute 'rtime'
sweep_number = 3
display(vol[sweep_number])
<xarray.Dataset>
Dimensions:                              (azimuth: 360, range: 324)
Coordinates:
    time                                 (azimuth) datetime64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
  * range                                (range) float32 250.0 ... 1.617e+05
  * azimuth                              (azimuth) float32 0.5274 1.53 ... 359.5
    elevation                            (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    latitude                             float32 ...
    longitude                            float32 ...
    altitude                             float32 ...
Data variables: (12/14)
    reflectivity                         (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    signal_to_noise_ratio                (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    reflectivity_vv                      (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    differential_reflectivity            (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    uncorrected_cross_correlation_ratio  (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    uncorrected_differential_phase       (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    ...                                   ...
    reflectivity_hh_clut                 (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    sweep_number                         int64 ...
    sweep_fixed_angle                    float32 ...
    sweep_mode                           |S32 ...
    pulse_width                          (azimuth) timedelta64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
    nyquist_velocity                     (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>

Clutter detection with Gabella

While in Switzerland, why not use the well-known clutter detection scheme by Marco Gabella et. al.

Wrap Gabella Clutter detection in Xarray apply_ufunc

The routine is implemented in wradlib in pure Numpy. Numpy based processing routines can be transformed to a first class Xarray citizen with the help of xr.apply_ufunc.

def extract_clutter(da, wsize=5, thrsnorain=0, tr1=6.0, n_p=6, tr2=1.3, rm_nans=False):
    return xr.apply_ufunc(
        wrl.clutter.filter_gabella,
        da,
        input_core_dims=[["azimuth", "range"]],
        output_core_dims=[["azimuth", "range"]],
        dask="parallelized",
        kwargs=dict(
            wsize=wsize,
            thrsnorain=thrsnorain,
            tr1=tr1,
            n_p=n_p,
            tr2=tr2,
            rm_nans=rm_nans,
        ),
    )

Calculate clutter map

Now we apply Gabella scheme and add the result to the Dataset.

swp = vol[sweep_number]
clmap = swp.reflectivity_hh_clut.pipe(
    extract_clutter, wsize=5, thrsnorain=0.0, tr1=21.0, n_p=23, tr2=1.3, rm_nans=False
)
swp = swp.assign({"CMAP": clmap})
display(swp)
<xarray.Dataset>
Dimensions:                              (azimuth: 360, range: 324)
Coordinates:
    time                                 (azimuth) datetime64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
  * range                                (range) float32 250.0 ... 1.617e+05
  * azimuth                              (azimuth) float32 0.5274 1.53 ... 359.5
    elevation                            (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    latitude                             float32 ...
    longitude                            float32 ...
    altitude                             float32 ...
Data variables: (12/15)
    reflectivity                         (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    signal_to_noise_ratio                (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    reflectivity_vv                      (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    differential_reflectivity            (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    uncorrected_cross_correlation_ratio  (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    uncorrected_differential_phase       (azimuth, range) float32 dask.array<chunksize=(360, 324), meta=np.ndarray>
    ...                                   ...
    sweep_number                         int64 ...
    sweep_fixed_angle                    float32 ...
    sweep_mode                           |S32 ...
    pulse_width                          (azimuth) timedelta64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
    nyquist_velocity                     (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    CMAP                                 (azimuth, range) bool dask.array<chunksize=(360, 324), meta=np.ndarray>

Plot Reflectivities, Clutter and Cluttermap

from osgeo import osr

wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
0
fig = plt.figure(figsize=(15, 12))

swpx = swp.sel(range=slice(0, 100000)).pipe(wrl.georef.georeference_dataset, proj=wgs84)

ax1 = fig.add_subplot(221)
swpx.reflectivity_hh_clut.plot(x="x", y="y", ax=ax1, vmin=0, vmax=60)
ax1.set_title("Reflectivity raw")

ax2 = fig.add_subplot(222)
swpx.CMAP.plot(x="x", y="y", ax=ax2)
ax2.set_title("Cluttermap")

ax3 = fig.add_subplot(223)
swpx.reflectivity_hh_clut.where(swpx.CMAP == 1).plot(
    x="x", y="y", ax=ax3, vmin=0, vmax=60
)
ax3.set_title("Clutter")

ax4 = fig.add_subplot(224)
swpx.reflectivity_hh_clut.where(swpx.CMAP < 1).plot(
    x="x", y="y", ax=ax4, vmin=0, vmax=60
)
ax4.set_title("Reflectivity clutter removed")
Text(0.5, 1.0, 'Reflectivity clutter removed')
../../_images/87d30faeaf35fc1e4b4d6e7ef8c8b5d24419cd4d71ca4ca07be7bc74e79d014e.png

SRTM based clutter and beamblockage processing

Download needed SRTM data

For the course we already provide the needed SRTM tiles. For normal operation you would need a NASA EARTHDATA account and a connected bearer token.

The data will be loaded using GDAL machinery and transformed into an Xarray DataArray.

extent = wrl.zonalstats.get_bbox(swpx.x.values, swpx.y.values)
extent
{'left': 7.545772342625747,
 'right': 10.120652596829073,
 'bottom': 45.14416090772705,
 'top': 46.9372157592286}
# apply fake token, data is already available
os.environ["WRADLIB_EARTHDATA_BEARER_TOKEN"] = ""
# set location of wradlib-data, where wradlib will search for any available data
os.environ["WRADLIB_DATA"] = "data/wradlib-data"
# get the tiles
dem = wrl.io.get_srtm(extent.values())
elevation = wrl.georef.read_gdal_values(dem)
coords = wrl.georef.read_gdal_coordinates(dem)
elev = xr.DataArray(
    data=elevation,
    dims=["y", "x"],
    coords={"lat": (["y", "x"], coords[..., 1]), "lon": (["y", "x"], coords[..., 0])},
)
elev
<xarray.DataArray (y: 2401, x: 4801)>
array([[ 422,  422,  422, ..., 1846, 1798, 1745],
       [ 422,  422,  422, ..., 1856, 1807, 1748],
       [ 422,  422,  422, ..., 1815, 1781, 1740],
       ...,
       [2403, 2379, 2357, ...,   14,   15,   15],
       [2408, 2387, 2365, ...,   14,   15,   14],
       [2428, 2411, 2397, ...,   13,   13,   13]], dtype=int16)
Coordinates:
    lat      (y, x) float64 47.0 47.0 47.0 47.0 47.0 ... 45.0 45.0 45.0 45.0
    lon      (y, x) float64 7.0 7.001 7.002 7.003 7.003 ... 11.0 11.0 11.0 11.0
Dimensions without coordinates: y, x

Plot Clutter on DEM

fig = plt.figure(figsize=(13, 10))
ax1 = fig.add_subplot(111)

swpx.CMAP.where(swpx.CMAP == 1).plot(
    x="x", y="y", ax=ax1, vmin=0, vmax=1, cmap="turbo", add_colorbar=False
)
ax1.set_title("Reflectivity corr")

ax1.plot(swpx.longitude.values, swpx.latitude.values, marker="*", c="r")

elev.plot(x="lon", y="lat", ax=ax1, zorder=-2, cmap="terrain")
<matplotlib.collections.QuadMesh at 0x7fb509825b50>
../../_images/21b5ed1d41c287b44da9f014dec957234c1153412fe7d967de518427f9cf1bc4.png

Use hvplot for interactive zooming and panning

Often it is desirable to quickly zoom and pan in the plots. Although matplotlib has that ability, it still is quite slow. Here hvplot, a holoviews based plotting framework, can be utilized. As frontend bokeh is used.

import hvplot
import hvplot.xarray

We need to rechunk the coordinates as hvplot needs chunked variables and coords.

cl = (
    swpx.CMAP.where(swpx.CMAP == 1)
    .chunk(chunks={})
    .hvplot.quadmesh(
        x="x", y="y", cmap="Reds", width=800, height=700, clim=(0, 1), alpha=0.6
    )
)
dm = elev.hvplot.quadmesh(
    x="lon",
    y="lat",
    cmap="terrain",
    width=800,
    height=700,
    clim=(0, 4000),
    rasterize=True,
)
dm * cl

Convert DEM to spherical coordinates

sitecoords = (swpx.longitude.values, swpx.latitude.values, swpx.altitude.values)
r = swpx.range.values
az = swpx.azimuth.values
bw = 0.8
beamradius = wrl.util.half_power_radius(r, bw)
rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(
    dem, nodata=-32768.0
)

rlimits = (extent["left"], extent["bottom"], extent["right"], extent["top"])
# Clip the region inside our bounding box
ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
rastercoords = rastercoords[ind[1] : ind[3], ind[0] : ind[2], ...]
rastervalues = rastervalues[ind[1] : ind[3], ind[0] : ind[2]]

polcoords = np.dstack([swpx.x.values, swpx.y.values])
# Map rastervalues to polar grid points
polarvalues = wrl.ipol.cart_to_irregular_spline(
    rastercoords, rastervalues, polcoords, order=3, prefilter=False
)

Partial and Cumulative Beamblockage

PBB = wrl.qual.beam_block_frac(polarvalues, swpx.z.values, beamradius)
PBB = np.ma.masked_invalid(PBB)
print(PBB.shape)
(360, 200)
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/qual.py:127: RuntimeWarning: invalid value encountered in sqrt
  numer = (ya * np.sqrt(a**2 - y**2)) + (a * np.arcsin(ya)) + (np.pi * a / 2.0)
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/qual.py:127: RuntimeWarning: invalid value encountered in arcsin
  numer = (ya * np.sqrt(a**2 - y**2)) + (a * np.arcsin(ya)) + (np.pi * a / 2.0)
CBB = wrl.qual.cum_beam_block_frac(PBB)
print(CBB.shape)
(360, 200)

Plotting Beamblockage

# just a little helper function to style x and y axes of our maps
def annotate_map(ax, cm=None, title=""):
    xticks = ax.get_xticks()
    ticks = (xticks / 1000).astype(int)
    ax.set_xticks(xticks)
    ax.set_xticklabels(ticks)
    yticks = ax.get_yticks()
    ticks = (yticks / 1000).astype(int)
    ax.set_yticks(yticks)
    ax.set_yticklabels(ticks)
    ax.set_xlabel("Kilometers")
    ax.set_ylabel("Kilometers")
    if not cm is None:
        plt.colorbar(cm, ax=ax)
    if not title == "":
        ax.set_title(title)
    ax.grid()
alt = swpx.z.values
fig = plt.figure(figsize=(15, 12))

# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)

# azimuth angle
angle = 270

# Plot terrain (on ax1)
ax1, dem = wrl.vis.plot_ppi(
    polarvalues, ax=ax1, r=r, az=az, cmap=mpl.cm.terrain, vmin=0.0
)
ax1.plot(
    [0, np.sin(np.radians(angle)) * 1e5], [0, np.cos(np.radians(angle)) * 1e5], "r-"
)
ax1.plot(sitecoords[0], sitecoords[1], "ro")
annotate_map(ax1, dem, "Terrain within {0} km range".format(np.max(r / 1000.0) + 0.1))
ax1.set_xlim(-100000, 100000)
ax1.set_ylim(-100000, 100000)

# Plot CBB (on ax2)
ax2, cbb = wrl.vis.plot_ppi(CBB, ax=ax2, r=r, az=az, cmap=mpl.cm.PuRd, vmin=0, vmax=1)
annotate_map(ax2, cbb, "Beam-Blockage Fraction")
ax2.set_xlim(-100000, 100000)
ax2.set_ylim(-100000, 100000)

# Plot single ray terrain profile on ax3
(bc,) = ax3.plot(r / 1000.0, alt[angle, :], "-b", linewidth=3, label="Beam Center")
(b3db,) = ax3.plot(
    r / 1000.0,
    (alt[angle, :] + beamradius),
    ":b",
    linewidth=1.5,
    label="3 dB Beam width",
)
ax3.plot(r / 1000.0, (alt[angle, :] - beamradius), ":b")
ax3.fill_between(r / 1000.0, 0.0, polarvalues[angle, :], color="0.75")
ax3.set_xlim(0.0, np.max(r / 1000.0) + 0.1)
ax3.set_ylim(0.0, 3000)
ax3.set_xlabel("Range (km)")
ax3.set_ylabel("Altitude (m)")
ax3.grid()

axb = ax3.twinx()
(bbf,) = axb.plot(r / 1000.0, CBB[angle, :], "-g", label="BBF")
axb.spines["right"].set_color("g")
axb.tick_params(axis="y", colors="g")
axb.set_ylabel("Beam-blockage fraction", c="g")
axb.set_ylim(0.0, 1.0)
axb.set_xlim(0.0, np.max(r / 1000.0) + 0.1)


legend = ax3.legend(
    (bc, b3db, bbf),
    ("Beam Center", "3 dB Beam width", "BBF"),
    loc="upper left",
    fontsize=10,
)
../../_images/3535121094292039f69af024059f9f96219b8c3827d46a73ce6b783f047d0e74.png

Plotting Beamblockage on Curvelinear Grid

Here you get an better impression of the actual beam progression .

def height_formatter(x, pos):
    x = (x - 6370000) / 1000
    fmt_str = "{:g}".format(x)
    return fmt_str


def range_formatter(x, pos):
    x = x / 1000.0
    fmt_str = "{:g}".format(x)
    return fmt_str
fig = plt.figure(figsize=(15, 8))

cgax, caax, paax = wrl.vis.create_cg(fig=fig, rot=0, scale=1)

# azimuth angle
angle = 270

# fix grid_helper
er = 6370000
gh = cgax.get_grid_helper()
gh.grid_finder.grid_locator2._nbins = 80
gh.grid_finder.grid_locator2._steps = [1, 2, 4, 5, 10]

# calculate beam_height and arc_distance for ke=1
# means line of sight
bhe = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=1.0)
ade = wrl.georef.bin_distance(r, 0, sitecoords[2], re=er, ke=1.0)
nn0 = np.zeros_like(r)
# for nice plotting we assume earth_radius = 6370000 m
ecp = nn0 + er
# theta (arc_distance sector angle)
thetap = -np.degrees(ade / er) + 90.0

# zero degree elevation with standard refraction
bh0 = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er)

# plot (ecp is earth surface normal null)
(bes,) = paax.plot(thetap, ecp, "-k", linewidth=3, label="Earth Surface NN")
(bc,) = paax.plot(thetap, ecp + alt[angle, :], "-b", linewidth=3, label="Beam Center")
(bc0r,) = paax.plot(thetap, ecp + bh0, "-g", label="0 deg Refraction")
(bc0n,) = paax.plot(thetap, ecp + bhe, "-r", label="0 deg line of sight")
(b3db,) = paax.plot(
    thetap, ecp + alt[angle, :] + beamradius, ":b", label="+3 dB Beam width"
)
paax.plot(thetap, ecp + alt[angle, :] - beamradius, ":b", label="-3 dB Beam width")

# orography
paax.fill_between(thetap, ecp, ecp + polarvalues[angle, :], color="0.75")

# shape axes
cgax.set_xlim(0, np.max(ade))
cgax.set_ylim([ecp.min() - 1000, ecp.max() + 2500])
caax.grid(True, axis="x")
cgax.grid(True, axis="y")
cgax.axis["top"].toggle(all=False)
caax.yaxis.set_major_locator(
    mpl.ticker.MaxNLocator(steps=[1, 2, 4, 5, 10], nbins=20, prune="both")
)
caax.xaxis.set_major_locator(mpl.ticker.MaxNLocator())
caax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(height_formatter))
caax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(range_formatter))

caax.set_xlabel("Range (km)")
caax.set_ylabel("Altitude (km)")

legend = paax.legend(
    (bes, bc0n, bc0r, bc, b3db),
    (
        "Earth Surface NN",
        "0 deg line of sight",
        "0 deg std refraction",
        "Beam Center",
        "3 dB Beam width",
    ),
    loc="lower left",
    fontsize=10,
)
../../_images/8d4cb06126c94b65fb4b16877c62c48c399e8e96ffbdf1192ea3bad62814e8db.png

Use Clutter and Beamblockage as Quality Index

Simple masking with cumulative beam blockage and Gabella.

swpx = swpx.assign({"CBB": (["azimuth", "range"], CBB)})
# recalculate georeferencing for AEQD
swpx = swpx.pipe(wrl.georef.georeference_dataset)
fig = plt.figure(figsize=(12, 4))

ax1 = fig.add_subplot(121)
swpx.reflectivity.plot(x="x", y="y", ax=ax1, cmap="turbo", vmin=0, vmax=60)
ax1.set_title(f"Signal Processor - {swpx.time.values.astype('M8[s]')}")
ax1.set_aspect("equal")

ax2 = fig.add_subplot(122)
# CBB > 0.5, CMAP == 1, RHOHV < 0.8 is masked
swpx.where(
    (swpx.CBB <= 0.5)
    & (swpx.CMAP < 1.0)
    & (swpx.uncorrected_cross_correlation_ratio >= 0.8)
).reflectivity_hh_clut.plot(x="x", y="y", ax=ax2, cmap="turbo", vmin=0, vmax=60)
ax2.set_title(f"Gabella+CBB+RHOHV - {swpx.time.values.astype('M8[s]')}")
ax2.set_aspect("equal")
../../_images/259220b32cf6af3c32b5f1ca3fe47be65d7756dfd6b53f3c7b993255b5eed8b5.png

Summary

We’ve just learned how to use \(\omega radlib\)’s Gabella clutter detection for single sweeps. Wrapping numpy based functions for use with xarray.apply_ufunc has been shown. We’ve looked into digital elevation maps and beam blockage calculations.

What’s next?

In the next notebook we dive into processing of differential phase.